function A = AutocorrelationInTreesIV(maxLag, colonySet, minPathLength)
%% AutocorrelationInTreesIV
%
% Author: Enrique Balleza
%% External params
%maxLag:                maximum lag to calculate
%colonySet:             colony structure or cell array with colony structures
%minPathLength:         only process path with lengths longer or equal to minPathLength
%
%
%% Internal params
% shrink:               number of points to delete in the lineages in order to reduce 
%                       experimental noise (these data are next to the gutters) 
%
%
%% Output
% A.Autocorr            <<Autocorr(tau)>> t and col average order doesn't matter.
%                       I take all the ~NaN matrix elements, add them together and 
%                       divide by the total number of ~NaN entries
%
% A.Autocorr_t_col      <<Autocorr(tau)>_t>_col
%
% A.colAutocorr         Matrix(nCol, tau) with rows the <<Autocorr(tau)>> of every colony
%                       in colonySet
%
% A.colAutocorrMedian   Median of A.colAutocorr
%
% A.colAutocorrMean     Mean of A.colAutocorr
%
% A.nData_per_col       Matrix(nCol, tau) in which rows are the number of data points used to
%                       calculate <<Autocorr(tau)>> of every colony. This
%                       matrix gives a relative contribution of the
%                       different colonies in the calculation of <<Autocorr(tau)>>
%% To do:
% Extend the selection of LoadColonies to include regexp
% Add the option to normalize by length 
% Add the other normalizations done by Paulsson

%% Decide between colony and colony set
%.......................................................
[colonySet, nCell, totalTime, maxTotalTime] = ColonyToColonySet(colonySet);

maxTotalTime


%Create set with all lineages from the different colonies
%.......................................................
[pathCellIndx, nLineage] = GetPathCellIndxSet(colonySet, minPathLength);


%Copy cell info to an array structure. This is translated into a faster code
%........................................................
nCol = length(colonySet);
cell_data = zeros(sum(nCell), maxTotalTime);
cell_dataL = zeros(1, length(cell_data(:,1)));
for i=1:nCol
    for j=1:nCell(i)
        indx = j+sum(nCell(1:i-1));
        cell_dataL(indx) = length(colonySet{i}.cell(j).flscnc.df);
        l = colonySet{i}.cell(j).size.l;
        df = colonySet{i}.cell(j).flscnc.df;
%         df_pos = df>=0;
%         df_neg = df<0;
%         df(df_pos) = sqrt(df(df_pos));
%         df(df_neg) = -sqrt(-df(df_neg));
        cell_data(indx,1:cell_dataL(indx)) = df;%./l;
    end
end


%Concatenate data for all different cell lineagesx
%.........................................................
lineage_data = NaN*zeros(sum(nLineage), maxTotalTime);
lineage_cellIndx = NaN*zeros(sum(nLineage), maxTotalTime);
lineage_lastIndx = NaN*zeros(sum(nLineage));
for i=1:length(pathCellIndx)
    
    for j=1:length(pathCellIndx{i})
        
        g_indxLineage = j + sum(nLineage(1:i-1));
        indxCell  = pathCellIndx{i}(j).lista;
        lastT = 0;
        for k=1:length(indxCell)
            g_indxCell = indxCell(k) + sum(nCell(1:i-1));
            n = cell_dataL(g_indxCell);
            dfT = cell_data(g_indxCell,1:n);
            lineage_data(g_indxLineage, lastT+1:n+lastT) = dfT;
            lineage_cellIndx(g_indxLineage, lastT+1:n+lastT) = indxCell(k);
            lastT = n+lastT;
        end
        lineage_lastIndx(g_indxLineage) = lastT;

    end
end


%Main Loop
%............................................
tic;
shrink = 3;

A = [];
A.Autocorr = zeros(1,maxLag);
A.Autocorr_t_col = zeros(1,maxLag);
A.colAutocorr = zeros(nCol,maxLag);
A.nData_per_col = zeros(nCol,maxLag);
A.colAutocorrMedian = zeros(1,maxLag);
A.colAutocorrMean = zeros(1,maxLag);

for lag=0:maxLag-1

    if mod(lag,10)==0, disp(['Estimating autocorrelation at lag: ' num2str(lag)]); end

    %Create contingency matrix
    %..................................................
    nTime = maxTotalTime-1-lag-shrink;         %calculate columns according to lag variable. Added -1 becasue last df point is always NaN
    tableData = NaN*ones(sum(nLineage), nTime);

    for col=1:nCol

        inf = 1+sum(nLineage(1:col-1));
        sup = sum(nLineage(1:col));
        rangeCol = inf:sup;        
        range2 = 1:(totalTime(col)-1-lag-shrink);

        tableData(rangeCol, range2) = ...
            AutocorrelationInTrees_Core(...
            pathCellIndx{col}, ...
            lineage_data(inf:sup, 1:totalTime(col)), ...
            lineage_cellIndx(inf:sup, 1:totalTime(col)), ...
            lineage_lastIndx(inf:sup), ...
            totalTime(col), nCell(col), lag, shrink);

    end

    A.Autocorr(lag+1) = sum(tableData(~isnan(tableData)))/( sum(sum(~isnan(tableData))) );

    for col=1:nCol
        inf = 1+sum(nLineage(1:col-1));
        sup = sum(nLineage(1:col));
        rangeCol = inf:sup;
        tableDataT = tableData(rangeCol,:);
        A.nData_per_col(col,lag+1) = sum(sum(~isnan(tableDataT)));
        A.colAutocorr(col,lag+1) = sum(tableDataT(~isnan(tableDataT)))/A.nData_per_col(col,lag+1);
    end

    AT = NaN*ones(1,nTime);
    nData = NaN*ones(1,nTime);
    for i=1:nTime
        here = ~isnan(tableData(:,i));
        nData(i) = sum(uint8(here));
        AT(i) = sum(tableData(here,i))/nData(i);
    end
    A.Autocorr_t_col(lag+1) = sum(AT)/length(AT); %unbiased estimate

end
toc;

for i=1:maxLag
    A.colAutocorrMedian(i) = median(A.colAutocorr(~isnan(A.colAutocorr(:,i)),i));
    A.colAutocorrMean(i) = mean(A.colAutocorr(~isnan(A.colAutocorr(:,i)),i));
end

end

%......................................................................
function [colonySet, nCell, totalTime, maxTotalTime] = ColonyToColonySet(colonySet)

if ~iscell(colonySet)
    colonySetT{1} = [];
    colonySetT{1} = colonySet;
    colonySet = colonySetT;
end
nCol = length(colonySet);
maxTotalTime = 0;
totalTime = zeros(1,nCol);

nCell = zeros(1,nCol);
for i=1:nCol
    totalTime(i) = length(colonySet{i}.bounds);
    nCell(i) = length(colonySet{i}.cell);
    if totalTime(i)> maxTotalTime, maxTotalTime = totalTime(i); end
end

end

%....................................................................
function [pathCellIndx, nLineage] = GetPathCellIndxSet(colonySet, minPathLength)

nCol = length(colonySet);
pathCellIndx{1} = [];
for i=1:nCol
    disp(['Creating redundant paths for colony ' num2str(i)]);

    pathCellIndx{i} = GetRedundantPathsII(colonySet{i});
    if minPathLength ~= 0
        pathCellIndxT = [];
        count = 1;
        for j=1:length(pathCellIndx{i})
            if length(pathCellIndx{i}(j).lista) >= minPathLength
               pathCellIndxT(count).lista =  pathCellIndx{i}(j).lista;
               count = count+1;
            end
        end
        pathCellIndx{i} = pathCellIndxT;
    end

end
nLineage = zeros(1,nCol);
for i=1:nCol, nLineage(i) = length(pathCellIndx{i}); end

end
